% skeleton experiment, what affects number of nodes in skeleton
function main()
MM=[500]; % 500
nn=[400]; % 300
khs=[1 5:5:45];
dr=[0.005:0.05:1.5005 1.755:0.25:3.005]; % percent extra, R=R*(1+dr)
objects=[1,2,3,4];
myss1=cell(length(objects),1);
myss2=cell(length(khs),1);

% change objects, change radius, static k,M,n,
for mm=1:length(objects)
	ss=[];
	for ii=1:length(MM)
		jj=1;kk=1;
		for ll=1:length(dr)
			ss1=exskel(MM(ii),nn(jj),khs(kk),dr(ll),objects(mm));
			ss=[ss ss1];
		end
	end
	myss1{mm}=ss;
end

% static object, change radius, change k, static M,n
for kk=1:length(khs)
    mm=1;
	ss=[];
	for ii=1:length(MM)
		jj=1;
		for ll=1:length(dr)
			ss1=exskel(MM(ii),nn(jj),khs(kk),dr(ll),objects(mm));
			ss=[ss ss1];
		end
	end
	myss2{kk}=ss;
end

figure, hold on
mycolors={'r','g','b','k','y'};
hs=[];
for ii=1:length(objects)
	hh=plot((1+dr)*100,myss1{ii},mycolors{ii});
	hs=[hs hh];
	fid = fopen(['o' num2str(objects(ii)) '.txt'],'wt');
	C=[(1+dr)*100;myss1{ii}];
	fprintf(fid,'%12.8f  %12.8f\n',C);
	fclose(fid);
end
legend(hs,'star 5','star 10','moon','ellips 3,1')
xlabel('radius % of smallest R')
ylabel('nodes in skeleton, original 300')
axis([100 400 0 300])


figure, hold on
hs=[];
for ii=1:length(khs)
	hh=plot((1+dr)*100,myss2{ii},'k-');
	hs=[hs hh];
	fid = fopen(['k' num2str(khs(ii)) '.txt'],'wt');
	C=[(1+dr)*100;myss2{ii}];
	fprintf(fid,'%12.8f  %12.8f\n',C);
	fclose(fid);
end
%legend(hs,'star 5','star 10','moon','circle','ellips 3,1')
xlabel('radius % of smallest R')
ylabel('nodes in skeleton, original 300')
axis([100 400 0 350])

return

function ss = exskel(M,n,kh,dr,otype)
nobj=1;
global nobj0;
global level;
global xctot;
global yctot;
global Rtot;
global isave;
isave = 0;

nobj0 = nobj;
centers = [[0,0]];
acc = 1e-10;


% Set the geometry.
contours = cell(nobj,1);
matsA    = cell(nobj,1);

R        = zeros(nobj,1);
xc       = zeros(nobj,1);
yc       = zeros(nobj,1);
ntot = n * nobj;

hmax = 0;
Ctot = [];
for i=1:nobj
	if (otype==1)
    	[C,len,xin,tmp] = get_geometry_star(n,centers(i,:).',5,0);
	elseif(otype==2)
		[C,len,xin,tmp] = get_geometry_star(n,centers(i,:).',10,0);
	elseif(otype==3)
		[C,len,xin,tmp] = get_geometry_moon(n,centers(i,:).',0);
	elseif(otype==4)
		[C,len,xin,tmp] = get_geometry_ellipsoid(n,centers(i,:).',3,1,0);
	end
    h_phy       = max(sqrt(C(2,:).^2 + C(5,:).^2));    
    hmax_phy    = max(hmax,h_phy);
    contours{i} = C;
    xc(i)       = mean(C(1,:));
    yc(i)       = mean(C(4,:));
    R(i)        = max(sqrt((C(1,:) - xc(i)).^2 + (C(4,:) - yc(i)).^2 ));
    Ctot        = [Ctot,C];          
end
xctot       = mean(Ctot(1,:));
yctot       = mean(Ctot(4,:));
Rtot        = max(sqrt((Ctot(1,:) - xctot).^2 + (Ctot(4,:) - yctot).^2 ));
% Create a reference problem
matA   = zeros(ntot);
% LOCAL_get_A_single_diag(Ctot,1:ntot,kh);
% Correct the diagonal blocks to take into account neighbors
for i=1:nobj
    for j=1:nobj
        i0 = (i-1)*n + 1; i1 = i0 + n - 1;
        j0 = (j-1)*n + 1; j1 = j0 + n - 1;
        if (i == j)
            matAii = get_A_single_diag(contours{i},1:n,kh);
            matsA{i}          = matAii;
        else
            matAii = get_A_offd_noquad(Ctot,i0:i1,j0:j1,kh);
        end
        matA(i0:i1,j0:j1) = matAii; 
    end
end
fprintf(1,'Object: %d   k: %f dr: %f\n',otype,kh,dr)
%[otype kh dr]
%[D,h_D]          = get_circle(contours,1,kh,2,0.1);
[D,h_D]          = get_circle2(contours,i,M,dr,kh);
[Is,Cs,V,matP,D] = get_skeleton(kh,contours{1},matsA{1},acc,D,h_D);
ss=size(Cs,2);
return
